Fish and health with national fish promotion programme

This is an open assessment looking at potential health effects of a national fish promotion program in Finland. The details of the assessment are described on Opasnet. This file contains the R code to run the assessment model.

Calculation is based on BONUS GOHERR project and its http://en.opasnet.org/w/Goherr_assessment.

# This is code Op_fi5923/model on page [[Kotimaisen kalan edistämisohjelma]]
# The code was forked from Op_fi5889/model on page [[Ruori]] and Op_en7748/model on page [[Goherr assessment]]
library(OpasnetUtils) # Install the newest version from https://github.com/jtuomist/OpasnetUtils not from CRAN.
library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# First empty all objects for a fresh start. Otherwise may be problems with CheckDecisions.
oempty(all=TRUE)
openv.setN(10)
dat <- opbase.data("Op_fi5923", subset="Malliparametrit")[-1]
dec <- opbase.data("Op_fi5923", subset="Decisions")[-1]
DecisionTableParser(dec)
CTable <- opbase.data("Op_fi5923",subset="CollapseMarginals")
for(i in 1:ncol(CTable)) {CTable[[i]] <- as.character(CTable[[i]])}
CollapseTableParser(CTable)
objects.latest("Op_en2031", code_name="subgrouping") # [[Exposure-response function]] subgrouping

## Loading objects:
##   subgrouping

cat("Laskennassa käytetty data.\n")

## Laskennassa käytetty data.

dat

cat("Tarkastellut päätökset.\n")

## Tarkastellut päätökset.

dec

cat("Aggregoidut marginaalit.\n")

## Aggregoidut marginaalit.

CTable

#' prepare adjusts the data table for ovariables. Requires function subgrouping from code Op_en2031/initiate on page [[Exposure-response function]]
#' @param dat data.frame
#' @param type type of data that is used. Must match content in column Type
#' @param drop columns to remove
#' @return data.frame
prepare <- function(dat, type=NULL, drop=NULL) {
  out <- dat
  if(!is.null(type)) out <- out[out$Type %in% type , ]
  if(!is.null(drop)) out <- out[!colnames(out) %in% drop]
  return(subgrouping(out))
}
amount <- Ovariable("amount", data = prepare(dat, "amount", c("Type","Response","Exposure_agent"))) # Filleted weight, i.e. no loss.
# Exposure:To child and To eater not needed, because dioxins are not (yet) included

conc_vit <- Ovariable(
  "conc_vit",
  ddata = "Op_en1838",
  subset = "Fineli data for common fish species"
)
#  dependencies = data.frame(Name="dummy"),
#  formula = function(...) {
#    df = opbase.data(
#      "Op_en1838", # [[Concentrations of beneficial nutrients in fish#Fineli]]
#      subset = "Fineli data for common fish species"
#    )
  df = conc_vit@data
  df$Nutrient[df$Nutrient=="D-vitamiini (µg)"] <- "Vitamin D"
  df$Nutrient[df$Nutrient=="rasvahapot n-3 moni-tyydyttymättömät (g)"] <- "Omega3"
  df$Nutrient[df$Nutrient=="rasvahappo 18:3 n-3 (alfalinoleenihappo) (mg)"] <- "ALA"
  df$Nutrient[df$Nutrient=="rasvahappo 22:6 n-3 (DHA) (mg)"] <- "DHA"
  df$Nutrient[df$Nutrient=="proteiini (g)"] <- "Fish"
  df$conc_vitResult[df$Nutrient=="Fish"] <- "1"
  df <- dropall(df[df$Nutrient %in% c("Vitamin D", "Omega3", "ALA", "DHA", "Fish") , ])
#  return(df)
#  }
#)
conc_vit@data <- df

exposure <- Ovariable(
  "exposure",
  dependencies = data.frame(Name=c("conc_vit", "amount")),
  formula = function(...) {
    
    # First, match KKE-classification with Fineli classification
    
    conc_vit$Kala <- NULL
    tmp <- Ovariable(
      output = data.frame(
        Kala = c("Kasvatettu", "Kaupallinen", "Kirjolohi", "Silakka", "Vapaa-ajan", "Muu tuonti", "Tuontikirjolohi", "Tuontilohi"),
        Fish = c("Whitefish", "Average fish","Rainbow trout", "Herring", "Average fish", "Average fish", "Rainbow trout", "Salmon"),
        Result = 1
      ))
      
    amount <- amount * 1000 / 5.52 /365.25 # M kg/a per 5.52M population --> g/d per person.
    out <- conc_vit * tmp * amount
    colnames(out@output)[colnames(out@output)=="Nutrient"] <- "Exposure_agent"
      
    return(out)
  },
  marginal = c(TRUE, TRUE, FALSE)
)
# Incidence-based data not needed
frexposed <- 1
population <- 0
incidence <- 0
ERFChoice <- Ovariable(
  "ERFchoice",
  data=data.frame(
    Response=c("Loss in child's IQ points","CHD2 mortality","Breast cancer","All-cause mortality","Depression"),
    Exposure_agent=c("DHA","Omega3","Omega3","Fish","Fish"),
    Result=1)
)
case_burden <- Ovariable("case_burden", data= prepare(dat,"case burden",c("Type","Exposure_agent","Unit")))
objects.latest("Op_en5917", code_name="InpBoD") # [[Disease risk]] InpBoD for GOHERR

## Loading objects:
##   InpBoD

InpBoD <- EvalOutput(InpBoD) # Evaluated because is not a dependency but an Input
InpBoD <- InpBoD[InpBoD$Country=="FI",colnames(InpBoD@output)!="Country"]
objects.latest("Op_en2261",code_name="BoDattr2") # [[Health impact assessment]]

## Loading objects:
##   BoDattr

tryCatch(BoDattr <- EvalOutput(BoDattr, verbose=TRUE))

##  Evaluating BoDattr ...
## 
##  - BoD fetched successfully!
## 
##  - PAF fetched successfully!
## - Evaluating BoD ...
## - - Evaluating case_burden ...
## 
##  done(0.02 secs)!
## - - Checking case_burden marginals ... Response, Iter, case_burdenSource recognized as marginal(s).

## Loading required package: reshape2

## - - Processing case_burden marginal collapses ... done!
## 
## - done(0.12 secs)!
## - Checking BoD marginals ... Response, Iter, BoDSource recognized as marginal(s).
## - Processing BoD inputs ... done!
## - Processing BoD marginal collapses ...

## Warning in oapply(variable, FUN = fun[[i]], cols = cols[[i]], na.rm = TRUE):
## While oapplying BoD, found NAs in indices: Group, InpBoDSource. They were
## automatically filled using fillna, which may result in a multiplied population.
## Please check your ovariable before using oapply.

##  done!
## - Evaluating PAF ...
## 
##  - - dose fetched successfully!
## 
##  - - ERF fetched successfully!
## 
##  - - RR fetched successfully!
## 
##  - - P_illness fetched successfully!
## 
##  - - sumExposcen fetched successfully!
## 
##  - - mc2d fetched successfully!
## 
##  - - mc2dparam fetched successfully!
## - - Evaluating dose ...
## 
##  - - - BW fetched successfully!
## - - - Evaluating exposure ...
## - - - - Evaluating conc_vit ...
## 
##  done(0 secs)!
## - - - - Checking conc_vit marginals ... Kala, Fish, Nutrient, Iter, conc_vitSource recognized as marginal(s).
## - - - - Processing conc_vit decisions ... done!
## - - - - Evaluating amount ...
## 
##  done(0 secs)!
## - - - - Checking amount marginals ... Kala, Scenario, amountSource recognized as marginal(s).
## 
## --- done(0.26 secs)!
## - - - Checking exposure marginals ... Fish, Exposure_agent, Iter, conc_vitSource, Adjust, Scenario, amountSource, exposureSource recognized as marginal(s).
## - - - Processing exposure marginal collapses ... done!
## - - - Evaluating BW ...
## 
##  done(0 secs)!
## - - - Checking BW marginals ... BWSource recognized as marginal(s).
## 
## -- done(25.94 secs)!
## - - Checking dose marginals ... Exposure_agent, Iter, conc_vitSource, Adjust, Scenario, amountSource, Scaling, exposureSource, BWSource, doseSource recognized as marginal(s).
## - - Processing dose marginal collapses ... done!
## - - Evaluating ERF ...
## 
##  - - - ERF_env fetched successfully!
## 
##  - - - ERF_omega3 fetched successfully!
## 
##  - - - ERF_mehg fetched successfully!
## 
##  - - - ERF_diox fetched successfully!
## 
##  - - - ERF_vit fetched successfully!
## 
##  - - - ERF_micr fetched successfully!
## 
##  - - - ERFchoice fetched successfully!
## - - - Evaluating ERF_env ...
## 
##  done(0.01 secs)!
## - - - Checking ERF_env marginals ... Exposure_agent, Response, Subgroup, Exposure, ER_function, Scaling, Exposure_unit, Observation, Iter, ERF_envSource recognized as marginal(s).
## - - - Evaluating ERF_omega3 ...
## 
##  done(0.01 secs)!
## - - - Checking ERF_omega3 marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_omega3Source recognized as marginal(s).
## - - - Evaluating ERF_mehg ...
## 
##  done(0 secs)!
## - - - Checking ERF_mehg marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_mehgSource recognized as marginal(s).
## - - - Evaluating ERF_diox ...
## 
##  done(0.01 secs)!
## - - - Checking ERF_diox marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_dioxSource recognized as marginal(s).
## - - - Evaluating ERF_vit ...
## 
##  done(0 secs)!
## - - - Checking ERF_vit marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, ERF_vitSource recognized as marginal(s).
## - - - Evaluating ERF_micr ...
## 
##  done(0 secs)!
## - - - Checking ERF_micr marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, ERF_micrSource recognized as marginal(s).
## - - - Evaluating ERFchoice ...
## - - - - Processing exposure marginal collapses ... done!
## 
## --- done(0.05 secs)!
## - - - Checking ERFchoice marginals ... Exposure_agent, ERFchoiceSource recognized as marginal(s).
## 
## -- done(2.86 mins)!
## - - Checking ERF marginals ... Exposure_agent, Response, Exposure, ER_function, Scaling, Observation, Iter, ERFSource recognized as marginal(s).
## - - Processing ERF marginal collapses ... done!
## - - Evaluating RR ...
## - - - Processing dose marginal collapses ... done!
## - - - Processing ERF marginal collapses ... done!
## 
## -- done(0.24 secs)!
## - - Checking RR marginals ... Exposure_agent, Response, ER_function, Scaling, Iter, ERFSource, conc_vitSource, Adjust, Scenario, amountSource, doseSource, RRSource recognized as marginal(s).
## - - Evaluating P_illness ...
## 
##  done(0 secs)!
## - - Checking P_illness marginals ... Response, Illness, Age, P_illnessSource recognized as marginal(s).
## 
## - done(6.1 mins)!
## - Checking PAF marginals ... Exposure_agent, Response, ER_function, Scaling, Iter, ERFSource, conc_vitSource, Adjust, Scenario, amountSource, doseSource, RRSource, PAFSource recognized as marginal(s).
## - Processing PAF marginal collapses ... done!
## 
##  done(6.93 mins)!
##  Checking BoDattr marginals ... Response, Iter, Group, InpBoDSource, Exposure_agent, conc_vitSource, Adjust, Scenario, amountSource, PAFSource, BoDattrSource recognized as marginal(s).

###################
# Graphs
trim <- function(ova) return(oapply(ova, NULL, mean, "Iter")@output)
plot_ly(trim(amount), x=~Scenario, y=~amountResult, color=~Kala, type="bar") %>%
  layout(yaxis=list(title="Kalan kokonaiskulutus Suomessa (milj kg /a)"), barmode="stack")

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

plot_ly(trim(conc_vit), x=~Nutrient, y=~conc_vitResult, color=~Kala, type="scatter", mode="markers") %>%
  layout(yaxis=list(title="Concentrations of nutrients (mg or ug /g)"))

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

tmp <- exposure / Ovariable(
  output = data.frame(
    Exposure_agent = c("Fish","Vitamin D", "Omega3", "ALA", "DHA"),
    Result = c(1, 1, 1000, 1000, 1000)
  ), marginal = c(TRUE, FALSE)
)
plot_ly(trim(tmp), x=~Scenario, y=~Result, color=~Exposure_agent, text=~Exposure_agent, type="bar") %>%
  layout(yaxis=list(title="Exposure to nutrients (g or ug /d)"))

summary(conc_vit, "mean")

summary(exposure,"mean")

summary(amount,"mean")

summary(BoD,"mean")

summary(BoDattr,"mean")

cat("Kalaperäisiä tautitaakkoja Suomessa\n")

## Kalaperäisiä tautitaakkoja Suomessa

tmp <- summary(oapply(BoDattr,NULL,sum,"Group")[BoDattr$Scenario=="BAU",])
data.frame(
  Altiste = tmp$Exposure_agent,
  Vaikutus = tmp$Response,
  Keskiarvo = as.character(signif(tmp$mean,2)),
  "95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
  Keskihajonta = signif(tmp$sd,2)
)#[rev(match(lev, tmp$Exposure_agent)),]

ggplot(trim(BoDattr), aes(x=Scenario, weight=BoDattrResult, fill=Response))+geom_bar()

ggplot(trim(BoDattr), aes(x=Scenario, weight=BoDattrResult, fill=Exposure_agent))+geom_bar()

plot_ly(trim(BoDattr), x=~Scenario, y=~BoDattrResult, color=~Response, text=~paste(Group, Exposure_agent, sep=": "), type="bar") %>%
  layout(yaxis=list(title="Disease burden (DALY /a); CHD2=coronary heart disease"), barmode="stack")

## Warning: Ignoring 12 observations

## Warning: n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning: n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

tmp <- summary(oapply(BoDattr,NULL,sum,c("Age","Exposure_agent")))
data.frame(
  Altiste = tmp$Response,
  Keskiarvo = signif(tmp$mean,2),
  "95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
  Keskihajonta = signif(tmp$sd,2)
)

################ Insight network
gr <- scrape(type="assessment")
objects.latest("Op_en3861", "makeGraph") # [[Insight network]]

## Loading objects:
##   makeGraph

gr <- makeGraph(gr)

## Loading required package: DiagrammeR

## Loading objects:
##   formatted
## Loading objects:
##   chooseGr

#export_graph(gr, "ruori.svg")
render_graph(gr)

##################### Diagnostics
objects.latest("Op_en6007", code_name="diagnostics")

## Loading objects:
##   showind
##   binoptest
##   showLoctable
##   ovashapetest

showLoctable()

showind()

## incidence is not an ovariable.
## population is not an ovariable.
## incidence is not an ovariable.
## frexposed is not an ovariable.
## sumExposcen is not an ovariable.
## mc2d is not an ovariable.
## mc2dparam is not an ovariable.